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Abstract 

In this article we provide a fast computational method in order to 
calculate the Moore-Penrose inverse of singular square matrices and 
of rectangular matrices. The proposed method proves to be much 
faster and has significantly better accuracy than the already proposed 
methods, while works for full and sparse matrices. 
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1 Introduction 

Let T be a n X n real matrix. It is known that when T is singular, then 
its unique generalized inverse T"'" (known as the Moore- Penrose inverse) is 
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defined. In the case when T is a real m x n matrix, Penrose showed that 
there is a unique matrix satisfying the four Penrose equations, called the 
generalized inverse of T. A lot of work concerning generalized inverses has 
been carried out, in finite and infinite dimension (e.g., [2l |8]). 

In a recent article [6], the first two authors provided a new method for the 
fast computation of the generalized inverse of full rank rectangular matrices 
and of square matrices with at least one zero row or column. In order to reach 
this goal, a special type of tensor product of two vectors was used, usually 
defined in infinite dimensional Hilbert spaces. In this work we extend our 
method so that it can be used in any kind of matrix, square or rectangular, 
full rank or not. The numerical experiments show that the proposed method 
is competitive in terms of accuracy and much faster than the commonly used 
methods, and can be also used for large sparse matrices. 
There are several methods for computing the Moore-Penrose inverse matrix 
(cf. [[2|). One of the most commonly used methods is the Singular Value 
Decomposition (SVD) method. This method is very accurate but also time- 
intensive since it requires a large amount of computational resources, espe- 
cially in the case of large matrices. On a recent work, Toutounian and Ataei 
^llj presented an algorithm based on the conjugate Gram-Schmidt process 
and the Moore-Penrose inverse of partitioned matrices, the CGS-MPi algo- 
rithm and they conclude that this algorithm is a robust and efficient tool 
for computing the Moore-Penrose inverse of large sparse and rank deficient 
matrices. 

In the present manuscript, we construct a very fast and reliable method 
(see the qrginv function in the Appendix) in order to estimate the Moore- 
Penrose inverse matrix. The computational effort required for the qrginv 
function (see Figure 1 and Tables 1, 4) in order to obtain the generalized 
inverse is substantially lower, particularly for large matrices, compared to 
those provided by the SVD and the method presented in ^Ij (CGS-MPi al- 
gorithm). In addition, we obtain reliable and very accurate approximations 
in each one of the tested cases (Tables 2, 3 and 4). In order to test this 
algorithm, we have used random singular matrices (see subsection 4.1) as 
well as a collection of singular test matrices (see subsection 4.2) with "large" 
condition number (ill-conditioned matrices) from the Matrix Computation 
Toolbox (see [9]). We also tested the proposed method on sparse matrices 
from the Matrix Market collection [7] (see subsection 4.3). In what follows, 
we make use of the high-level language Matlab both for calculations of the 
generalized inverse matrix, as well as for testing the reliability of the ob- 
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tained results. Specifically, the Matlab 7.4 (R2007a) Service Pack 3 version 
of the software was used on an Intel Core i7 920 Processor system running 
at 2.67GHz with 6 GB of RAM memory using the Windows XP Professional 
64-bit Operating System. 



2 Preliminaries and notation 

We shall denote by M*"^" the linear space of all m x n real matrices. For 
T G M™'^"', the generalized inverse G R*^^*^ (known as the Moore- Pen- 
rose inverse) is the unique matrix that satisfies the following four Penrose 
equations: 

where T* denotes the transpose matrix of T. The number r = dim7?.(T) is 
called the rank of T and ( , ) denotes the usual inner-product in M". 
According to [TU], for each x G M'^, we consider the mapping 

e (g) / : M'^ ^ M'^ with (e O f){x) = (x, e)f. 

Assume that {ei,...,en} and {/i,...,/n} are two collections of orthonormal 
vectors and linearly independent vectors of M.^ with n < k, respectively. 
Then, every rank-77, operator T can be written in the form T = X]r=i ^i® fi- 
We shall refer to this type of tensor product as the tensor-product of the 
collections {ei,...,e„} and {/i,.-.,/n}- The adjoint operator T* of T is the 
rank-n operator T* = XlILi fi® ^i- 

The tensor-product of two collections of vectors, as defined above, is a 
linear operator, therefore, it has a corresponding matrix representation T. 
We shall refer to this matrix T as the tensor-product matrix of the given col- 
lections. In order to compute the Moore-Penrose inverse of the corresponding 
tensor-product matrix, we use the following theorem. 

Theorem 2.1 ([5J, Theorem 3.2) LetH be a Hilbert space. IfT = Er=i 
fi is a rank-n operator then its generalized inverse is also a rank-n operator 
and for each x El-L, it is defined by the relation 

n 

T^x = ^ Xi{x)ei, 

i=l 

where the functions Aj are the solution of an appropriately defined nxn linear 
system. 
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3 The computational method 



In [6], based on theorem 1, the authors developed an algorithm (the ginv 
function) for computing the generalized inverse of full rank matrices, and of 
square matrices with at least one zero row or column and the rest of the 
matrix full rank. In other words, our main concern was to calculate the 
corresponding Aj in the expansion 

n 

T^x = ^ Xi{x)ei, 

i=l 

where {ei, e„} are the first n vectors of the standard basis of M'^, in order 
to provide the generalized inverse 

To extend this result for any kind of matrix, we will make use of the QR 
factorization, as well as the reverse order law for generalized inverses. The 
following proposition is a restatement of a part of R. Bouldin's theorem [1] 
which holds for operators and matrices (see also [3], 

Proposition 3.1 Let A, B be bounded operators on with closed range. 
Then {ABy = if and only if the following three conditions hold: 

i) The range of AB is closed, 

a) A'^A commutes with BB* , 

Hi) BB^^ commutes with A* A. 

Using the above proposition, we have the following result, which can be found 
also, in a similar form but with a different proof, in |2]. 

Proposition 3.2 Let A = QR be the QR factorization of A. Then, ^4^ = 
R^Q*. 

Proof: We must prove that the conditions of Bouldin's theorem hold. The 
first condition always holds, since in the case of matrices the range is always 
closed. For the second condition, it is easy to see that since Q is a unitary 
matrix, Q"^ = Q* = Q~^ and so 

Q^QRR* = Q-^QRR* = IRR* = RR*I = RR*Q^Q. 

The third condition can be proved in a similar way. 

Using the QR factorization, the matrix R is upper triangular but not nec- 
essarily of full rank. So a variant of the QR method must be used, the QR 
with column pivoting as described in the following form from [T^ : 
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Theorem 3.3 ([12], Theorem 3.3.11) Let A G M"'''" matrix, withrank{A) 
r > 0. Then there exist matrices A,Q,R such that A is obtained by A by 

Rn Ri2 


is nonsingular and upper triangular and A = QR. 



permuting its columns, Q G 
Rn G W'' 



is orthogonal, R 



Using the above theorem, we have that AP = QR, where P is a per- 
mutation matrix (therefore unitary). By proposition 13.21 we have that = 
PR^Q*. 

To calculate the rank of Ru, one needs only the number of its columns 
that have at least one value above a tolerance level in absolute terms. This 
tolerance is set equal to 10~^, which is also used by Toutounian and Ataei 
[11], and turns out to provide accurate results. The implementation of all 
the above ideas are presented in the qrginv function (see the Appendix). 



4 Numerical experiments 

In this section we perform numerical experiments comparing Matlab's pinv 
function, Toutounian and Ataei' s method [TT] ( CGS-MPi algorithm) and the 
proposed method qrginv function. Testing of pinv, CGS-MPi and qrginv 
was performed separately for random singular matrices and for singular test 
matrices with "large" condition number from the Matrix Computation Tool- 
box (see [9]). We also tested the proposed method in sparse matrices and we 
obtained very fast and accurate results. 



4.1 Random singular matrices 

We are comparing the performance of the proposed method (qrginv) to that 
of the SVD method used by Matlab (pinv) and the method proposed in 
[TT] (CGS-MPi algorithm), on a series of random singular matrices with rank 
2", n = 8, 9, 10, 11, 12. In addition, the accuracy of the results is examined 
with the matrix 2-norm in error matrices corresponding to the four properties 
characterizing the Moore-Penrose inverse. 

In Table 1 we can see the time efficiency for the same matrices of the 
proposed method, the CGS-MPi algorithm and the pinv command used by 
Matlab. The qrginv method needs about 4% up to 11% the corresponding 
time needed by the pinv function to calculate the Moore-Penrose inverse, 
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Table 1: Computational Time Results; Random singular matrices 



Rank 


pinv 


CGS-MPi 


qrginv 


28 


0.496 


1.281 


0.0247 


29 


1.185 


15.183 


0.104 


210 


14.066 


448.172 


1.556 


2" 


152.384 


11374.47 


9.259 


2I2 


1081.81 




45.061 



Notes: Time is measured in seconds. CGS-MPi 
algorithm was not able to produce numerical 
results for a rank 2^^ matrix, even after 3 days. 

depending on the matrix. On the other hand the CGS-MPi turns to be com- 
putantionally demanding requiring from 50 times up to more than 1200 times 
the corresponding time needed by qrginv. Furthermore, the larger the rank 
of the matrix, the greater this difference, so that for a matrix with rank 2^^, 
the CGS-MPi algorithm was not able to produce a numerical result, even after 
3 days running, while qrginv needs only up to 45 seconds. 

In Table 2, the accuracy of the proposed method is tested based on the 
2-norm errors. It is evident that the proposed method (qrginv) produced a 
reliable approximation in all the tests that were conducted. The associated 
errors are greatly lower than the corresponding of the other methods, while in 
certain (many) cases are equal to zero, under the 64-bit IEEE double precision 
arithmetic in Matlab. Therefore, the proposed method allows us for a both 
fast and accurate computation of the Moore-Penrose inverse matrix. 

4.2 Singular test matrices 

In this section we use a set of singular test matrices that includes 13 singular 
matrices, of size 200 x 200, obtained from the function matrix in the Matrix 
Computation Toolbox [9j (which includes test matrices from Matlab itself). 
The condition number of these matrices range from order 10^^ to 10^^^. For 
comparative purpose we also apply as in the previous section, Matlab's pinv 
function which implements the SVD method and CGS-MPi algorithm. Since 
these matrices are of relatively small size and so as to measure the time 
needed for each algorithm to compute the Moore-Penrose inverse accurately, 
each algorithm runs 100 distinct times. The reported time is the mean time 
over these 100 replications. The error results are presented in Table 3, while 
the time responses are shown in Figure 1. 



6 



Table 2: Error Results; Random singular matrices 

Method Rank ||TTtr-r||2 ||Ttrrt - rt||2 \\TT^ - {TT^)*\\2 ||Ttr - (Ttr)* ||2 



pinv 1.21 X 10-12 5.58 x 10"" 2.98 x 10"^^ 3.46 x lO"^'^ 

CGS-MPi 2* 2.13 x 10"* 7.44 x lO"'' 2.01 x 10"" 5.29 x 10"'^ 



qrginv 



pmv 



qrgmv 



1.44 X 10-13 



3.07 x 10-12 2.31 X 10-12 1.49 x IQ-i^ 7.33 x 10" 



CGS-MPi 29 1.62 x 10-^ 9.79 x 10-^ 9.68 x IQ-" 5.18 x 10-^ 



3.62 x 10-13 



pinv 8.24 x lO-i^ 3.18 x lO-i^ 1.63 x IQ-i^ 1.62 x lO-i^ 

CGS-MPi 2io 1.57 x 10-^ 3.70 x lO-'' 2.90 x 10-1° 1.65 x 10-'^ 

qrginv 7.47 x 10-" 



pinv 4.11 x 10-11 2.90 X 10-" 9.32 x IQ-i^ 9.45 x lO-i^ 

CGS-MPi 2" 6.99 x lO-'^ 4.19 x 10-^^ 5.66 x 10-^ 1.15 x 10-^^ 

qrginv 6.31 X 10-12 



pinv 1.19 x 10-1° 7.52 x 10-1° 2.43 x IQ-n 5.69 x 10-" 

CGS-MPi 2i2 - 

qrginv 1.01 X 10-" 



Notes: Zeros (0) denote numbers close to zero under 64-bit IEEE double precision arithmetic. The CGS-MPi 
algorithm was not able to produce numerical results for matrix with rank 2i2 even after 3 days running. 



The results show a clear ordering of the three methods for this set of 
test problems, with qrginv in the first place, followed by pinv and then 
by CGS-MPi algorithm. The worst cases for pinv, with comparatively large 
errors for the \\T^TT^ ~ ^^^Ib norm, are the lotkin, prolate, hilb, and vand 
matrices. The CGS-MPi algorithm has very large errors in the cases of Kahan, 
lotkin, prolate, hilb, magic and vand matrices, for all tested norms. Moreover, 
the algorithm failed to produce a numerical result for the chow matrix, since 
the computed Moore-Penrose inverse included only NaNs. Nevertheless, we 
observe that there are also cases (in certain norms) that the CGS-MPi algo- 
rithm has smaller error than pinv . On the other hand, the proposed qrginv 
method gives very accurate results for all matrices and proves to be overally 
more efficient than the other two methods. 
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Figure 1: Time efficiency for the pinv, qrginv functions and the CGS-MPi 
algorithm 

4.3 Matrix Market sparse matrices 

In this section we test the proposed algorithm on sparse matrices, from the 
Matrix Market collection [7]. We follow the same method and the same 
matrices as in [TT], while the matrices are taken as rank deficient: A_Z = 
[A Z], where A is one of the chosen matrices, Z is a zero matrix of order 
m X 100 and m takes values shown in Table 4, as in 

The proposed algorithm is tested against the CGS-MPi algorithm, since 
this method is proposed by Toutounian and Ataei [TT] as suitable for large 
and sparse matrices. The pinv function of Matlab is not applicable in sparse 
matrices and thus ommited. We observe that in sparse matrices as well, the 
proposed method seems to greatly outperform the CGS-MPi algorithm, both 
in terms of speed and accuracy. 

5 Concluding Remarks 

We proposed a new method for calculating the Moore- Penrose inverse of sin- 
gular square, rectangular, full or sparse matrices. It is apparent that the 
proposed method provides a substantially faster numerical way for calculat- 
ing the Moore- Penrose inverse of a given matrix (see Figure 1 and Tables 1, 
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Table 3: Error Results; Singular test matrices 









II J J i- -i ||2 


IITtTTt 


— Tt lU 
112 


\\l 1 (J. J. ) ||2 




^ ) 112 


chow 


8.01849 X 10^^^ 


pinv 


4.741 X 10-^3 


1.896 


X 


10-13 


2.313 X 10-13 


2.261 X 


10-13 


(r=199) 




ULio-JViJ-'i 




















qrginv 



















cycol 


2.06 X lO"^ 


pinv 


1.539 X 10"^^ 


1.637 


X 


10- 1« 


4.285 X 10-1^ 


3.996 X 


10-15 


(r— 5U; 




OvjO-iVlJr'l 


'\ f^n7 V 1 n~i4 

O.OU/ X lU 


6.365 


X 


10-17 

i-U 


u.oyi A ±\j 


1.236 X 


10-15 






qrginv 











Q 







gearmat 


3.604 X lO^'' 


pinv 


1.899 X 10-" 


3.033 


X 


10-13 


8.063 X 10-14 


7.561 X 


10-14 


(r = 199) 






z.o ( y X lu 




X 




fi 1 1 V 1 n— 14 

D.llO X LU 


5.478 X 


10-13 






qrginv 


1 Q9^ y in~^^ 

















kahan 


2.30018 X 10^* 


pinv 


1.432 X 10-^3 












9.070 X 


10-14 


(r = 199) 






4 ni 1 V "1 n+i^ 

^.U-Li A lU 




X 






5.765 X 


10+18 






qrginv 


fi Qfi4 X 1 D"^'"' 

\J.iJ\J'-t. /\ _LU 

















lotkin 


8.97733 X 10^^ 


pinv 


0.0001 


lUo.o 


X 


10+6 


0.001 


0.0011 




(- 1 Q\ 




o o- ivi jr 1 


9Q V 1 n+'^ 

-OU A _LU 


1.252 


X 


10+14 


Q 


74.5 X 10+6 






qrginv 


O.O'iU A J_U 

















prolate 


5.61627 X 10^'' 




0.0002 


7.398 


X 


10^ 


0.0096 


0.0081 




(r = 117) 




CGS-MPi 


2.218 X 10+1'' 


7.885 


X 


10+30 


0.625 


4.796 X 


10+15 






qrginv 


6.588 X 10-" 

















hilb 


1.17164 X 10^^ 


pinv 


0.0001 


3.184 


X 


10" 


0.0033 


0.0051 




(r = 20) 




CGS-MPi 


39.6 X 10+5 


2.811 


X 


10+12 


5.533 X 10-1° 


30.95 X 


10+5 






qrginv 


6.941 X 10-12 

















magic 


4.92368 X 10^^ 


pinv 





2.491 


X 


10-19 


9.701 X 10-1* 


4.236 X 


10-14 


(r = 3) 




CGS-MPi 


1.755 X 10+18 


8.333 


X 


10+19 


293.29 


6.329 X 


10+14 






qrginv 


4.566 X 10-1" 

















vand 


1.16262 X 10^" 


pinv 


0.0003 


179.8 


X 


10+« 


0.0022 


0.0019 




(r = 34) 




CGS-MPi 


6.498 X 10+11 


5.299 


X 


10+20 





2.39 X 10+11 






qrginv 


5.553 X 10-11 


















Notes; Zeros (0) denote numbers close to zero under 64-bit IEEE double precision arithmetic. In 
parenthesis is denoted the rank (r) of each matrix. The CGS-MPi algorithm was not able to produce 
numerical results for the chow matrix. 



4). Also, it is evident, from Tables 2, 3 and 4 that the proposed function 
(qrginv)provides a far more reliable approximation in all the tests, compared 
to other existing methods. 

Appendix: Matlab code of the qrginv function 

function qrginv = qrginv(B) [N,M] = size(B); [Q,R,P] = qr(B); 
r=sum(any(abs(R)>le-5,2)) ; Rl = R(l:r,:); R2 = ginv(Rl); R3 = [R2 
zeros (M, N-r) ] ; A = P*R3*Q'; qrginv = A; 
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Table 4: Error and Computational Time Results; Matrix Market sparse ma- 
trices 



Matrix 


Time 


Method 


\\TT^T — T\\2 


\\T^TT^ — II2 


||TTt — (TTt)*||2 


IIT+T — (T+T)* II2 


WELL1033_Z 


0.0716 


CGS-MPi 


2.801 X 10"" 





9.307 X 10-12 


1.303 X 10-1° 


(m = 1033) 


0.0111 


qrginv 


1.142 X 10-12 











WELL1860.Z 


0.4626 


CGS-MPi 


2.184 X 10-11 


1.136 X 10-1° 


7.135 X 10-12 


7.236 X 10-11 


(m = 1860) 


0.0315 


qrginv 


1.89 X 10-12 











ILCC1850_Z 


0.4546 


CGS-MPi 


61.106 X 10+4 


1.687 X 10+8 


2.637 X 10-1° 


50.971 X 10+4 


(m = 1850) 


0.0315 


qrginv 


5.809 X 10-11 











GR-30-30.Z 


0.9618 


CGS-MPi 


4.765 X 10-1° 


5.91 X 10-11 


2.206 X 10-11 


6.746 X 10-1° 


(m = 900) 


0.0336 


qrginv 


7.321 X 10-12 











WATTl.Z 


1.2687 


CGS-MPi 


8 


101.6 X 10+« 


7.743 X 10-i« 


27.783 


(m = 1856) 


0.0064 


qrginv 















Notes; Zeros (0) denote numbers close to zero under 64-bit IEEE double precision arithmetic. In 



parenthesis is denoted the row size of each matrix (m x 100) . Time is measured in seconds. 

In case the matrix of interest is sparse, the third line of the code is replaced 
with 

[Q,R,P] = spqr(B); 

since the function qr embeded in Matlab does not support sparse matrices 
under this format. The function spqr is part of the SuiteSparse toolbox, built 
by Professor Timothy A. Davis, University of Florida and can be downloaded 
electronically from hUp://www. cise.ufl. edu/research/sparse/SuiteSparse/. 
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